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ABSTRACT 

In a previous work Lazarian and Pogosyan suggested a technique to extract velocity and den- 
sity statistics, of interstellar turbulence, by means of analysing statistics of spectral line data 
cubes. In this paper we test that technique, by studying the effect of correlation between ve- 
locity and density fields, providing a systematic analysis of the uncertainties arising from the 
numerics, and exploring the effect of a linear shear We make use of both compressible MHD 
simulations and synthetic data to emulate spectroscopic observations and test the technique. 
With the same synthetic spectroscopic data, we also studied anisotropics of the two point 
statistics and related those anisotropics with the magnetic field direction. This presents a new 
technique for magnetic field studies. The results show that the velocity and density spectral 
indices measured are consistent with the analytical predictions. We identified the dominant 
source of error with the limited number of data points along a given line of sight. We decrease 
this type of noise by increasing the number of points and by introducing Gaussian smoothing. 
We argue that in real observations the number of emitting elements is essentially infinite and 
that source of noise vanishes. 

Key words: turbulence - ISM: general, structure - MHD - radio lines: MHD. 



1 INTRODUCTION 

The very large Reynolds numbers (defined as the ratio of the in- 
ertial force to the viscous force acting on a parcel of gas) clearly 
suggest that the interstellar medium is turbulent. Understanding of 
this turbulence is crucial to the correct and complete description 
of many physical processes that take place in the ISM: molecular 
cloud dynamics, star formation, heat transfer, magnetic reconnec- 
tion, accretion disks, cosmic ray propagation and diffusion, just to 
mention a few. 

To adequately describe turbulence we need to use statistical 
methods, to extract the underlying regularities and reject incidental 
details. So far, studies of turbulence statistics have been most suc- 
cessful using interstellar scintillations, and have provided impor- 
tant information of the density statistics on scales 10*^ — 10^ cm, 
(see iNaravan & GoodmanllT989l: ISoangler & Gwinj|l99(]l) . These 
studies are done with the ionised media, and are restricted to the 
density fluctuations. We expect to observe power-law power spec- 
trum of density in a turbulent medium, however density fluctu- 
ations alone fail to distinguish between fossil and active turbu- 
lence. The distribution of sizes of sand grains in a beach follows 
a power-law but it is not turbulence. Therefore would be better to 
study a more direct and dynamically important statistic: velocity 
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fluctuations. Studies of the neutral media have included velocity 
information, e.g. using spectral line widths, centroids of veloc- 
itv (Mies ch & Bally.. 1994' .Scalo.. 1 987.): or more fancy statistics 
faever & Schloerblll997t iRosolowskv et alJll999l: fsrunt & HeveJ 
l2002jd ). 

Turbulence statistics is a very important link between the- 
ory and observations. Recent breakthroughs in turbulence studies 
have been possible with the use of scaling laws (see reviews by 
Cho & Lazarian 2003: Cho, Lazarian & Vishniac 2003, and refer- 
ences therein). Scaling laws predictions can be confronted with ob- 
servations. 

A discussion of various approac hes to the study of turbulence 
with spectral line data can be found in Lazarian ( 1999). In particular 
the problem of the contribution of velocity a nd density fluctuations 
to the emissivity statistics was addressed in lLazarian & Poeosvanl 
hereafter LPOO), where analytical description of the power 
spectrum of the emissivity in velocity channels was obtained. This 
study provides a technique to extract statistics of density and ve- 
locity from spectral line data cubes. H I has been chosen as a test 
case because it was possible to disregard self-absorption, and con- 
sider the emissivity linearly depending on the density (see review 
in Lazarian, Pogosvan & EsQuivel 2002). The technique, termed as 
velocity-channel analysis (or simply, VGA), is based on the vari- 
ations of the power spectra in velocity channels with the change 
of velocity resolution. Predictions in LPOO have been confirmed 
through observations iStanimirovic & LazariarJ r200lh . The first 
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tests using numerical data in lLazarian et alj i200lh also confimied 
the predictions in LPOO. However they did not address the problems 
of numerical noise, effects of magnetic field and galactic shear. 
Moreover the study of velocity-density correlations was performed 
over a very limited dynamical range. 

In this work we subject the predictions in LPOO to further 
scrutiny. In particular, we systematicall y analyse the effect s of 
velocity-density correlations, and unlike Lazarian et al. (2001) we 
include a study of the effect of shear. In addition we address the is- 
sue of an anisotropic turbulent cascade, with their implications on 
the VCA. We also analyse the sources of numerical noise. 

In §2 we review the basic problem, we then (§3) test the 
extraction of the spectral indices using the VCA, and obtain the 
density-velocity correlations. In §4 we study the effect of shear in 
the VCA. The anisotropic cascade, and a possible new method to 
extract information of the magnetic field direction can be found in 
§5. And we finally provide a summary of our results (§6). 



2 STATISTICS OF SPECTRAL LINE DATA AND VCA 

Spectral line observations contain important information of den- 
sity and velocity. The problem is that both velocity fluctuations and 
density fluctuations contribute to the intensity of a spectral line at a 
given velocity, and their separation is far from trivial. 



2.1 3D turbulence statistics 

Turbulence, due to its stochastic nature is studied using statis- 
tical tools, such as 3D correlation and structure functions (see 
iMonin & Yaglomll975h . For instance, the density correlation func- 
tion can be expressed as 



C(r) = {p{x)p{x + r)), 



(1) 



where x is the spatial position (xyz in Cartesian coordinates), r 
is a spatial separation (or 'lag'), and the angular brackets denote 
an average to be performed over all x space. The power spectrum 
provides an alternative description and is related to the correlation 
function (CF) as 



P{k) = / dre*'="'e(r) 



(2) 

where k is the wave-number, related to the scale under study like 
k — 27r/r, and the integration is performed over all 3D space. For 
power-law statistics, the A'^ dimensional power spectrum (Pm oc 
A:") and the correlation (^m oc r'") have indices related as n = 
— A'^ — m, in other words^ 

{spectral index) — {—spatial dimensions — CF index). (3) 
To describe velocity fields is customary to use structure functions: 



SF{r) = {[v{x) -v{x + r)Y). 



(4) 



They provide a more formal treatment of velocity fields 
(Monin & Yaalom 1975). In the case of isotropic velocity fields 
it differs from the correlation function only by a constant. There- 
fore for power-law statistics it has the same spectral index m as 



^ Notice that the spectral indices contain the sign. This notation is unfor- 
tunate (contrary to a conventional spectral index of the form P]v oc fc^", 
with n positive). However, it was introduced in LPOO, and for consistency 
we adhere to it. 



the correlation function (i. e. SF oc r™). For instance, Kolmgorov 
turbulence scales as n oc r^^^, which corresponds to m = 2/3. 



2.2 Statistics in position-position-velocity (PPV) space 

In spectroscopic observations we don't observe the gas distribution 
in real space coordinates x = {x, y, z). Instead the intensity of the 
emission in a given spectral line is measured towards some direc- 
tion in the sky, given by X = (x, y), and at a given line-of-sight 
(LOS) velocity v. Thus, the coordinates of observational data are 
{X, v). Therefore, observational data are contained in, so called, 
position-position- velocity (or simply 'PPV') cubes. If we identify 
the z coordinate with the LOS, then the relation between the real 
space and PPV descriptions is that of a map (X, z) {X ,v). 

At a given position x the LOS component of the velocity v can 
be decomposed as a regular flow velocity Ureg{x), plus a thermal 
component Vthermai and a turbulent component u{x). This way the 
observed Doppler shifted atoms follow a Maxwellian distribution 
of the form 



(l)v{x)dv ■ 



(27r/3)i/2 



exp 



[V - Vreg{x) - u{x)f 

2/3 



d^;, (5) 



where (3 = KsT/m, m is the atomic mass, T the gas temperature, 
and Kb the Maxwell-Boltzmann constant. 

The density of emitters in PPV space ps{X,v) can be ob- 
tained by integration along the LOS: 



Ps{X,v)dXdv = 



dz p{x)(l>{x) 



dXdv, 



(6) 



where p{x) is the mass density of the gas in spatial coordinates. 
This expression just counts the number of atoms along the LOS, at 
a given position X with a z component of velocity in the interval 
[v,v + dv], and the limits of integration are defined by the extent 
of spatial distribution of emitting gas. 

The statistics available in PPV are the correlation functions, 
= {ps{Xi,vi)ps{X2,V2)), and the power spectra (this can be 
ID, 2D or 3D) of emissivity fluctuations. Relations between this 
statistics and the underlying velocity and density spectra were es- 
tablished in LPOO, where the 2D power spectrum in velocity chan- 
nels was used. This spectrum can be directly obtained with radio- 
interferometric observations (see iLazariarh 1 995i) . 



2.3 Velocity-channel analysis (VCA) 

This technique was developed in LPOO to extract the velocity and 
density spectral indices of turbulence from spectral line data. LPOO 
provide analytical predictions for the emissivity power spectral in- 
dices in velocity channels, as a function of the velocity resolution 
employed. The relative contribution to the intensity in a velocity 
channel (also referred as 'slice') to the total intensity fluctuations 
changes in a regular fashion with the width of the velocity slice. 
This is easy to understand qualitatively since more velocity fluctua- 
tions within a channel are averaged out as we increase the thickness 
of the channel, decreasing the relative contribution of velocity. 

The notation that we will use is such that: n is the 3D density 
spectral index (Pp oc fc"); jj, the 3D velocity spectral index (Pv oc 
fc^); 7 corresponds to the spectral index in velocity slices (P„ oc 
fc^); and m is the structure function index of the velocity. From eq. 
jS) is clear that 



P - 



'3 — m. 



(7) 
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Table 1. A summary of analytical results derived in LPOO. 



Slice 


Shallow 3-D density Steep 3-D density 


thickness 


P„ (X fc", n > -3 Pn oc fe", n < -3 


2-D intensity spectrum for thin'* slice 




2-D intensity spectrum for thick'' slice 




2-D intensity spectrum for very thick'^ slice 





channel width < velocity dispersion at the scale under study 
channel width > velocity dispersion at the scale under study 
substantial part of the velocity profile is integrated over 



For a given 3D density power spectrum two distinct regimes 
are present when: a) n > —3 ('shallow' density case), and (b) 
n < —3 ('steep' density case). A summary of the predictions in 
LPOO for the intensity spectral index in 2D velocity channels is 
given m Table □ In both cases (shallow or steep) the power-law 
index 'gradually steepens' with the increase of the velocity slice 
thickness. In the thickest velocity channels all the velocity infor- 
mation is averaged out and we naturally get the density spectral 
index n. On the other hand, the velocity fluctuations dominate in 
thin channels. For instance, in the case of a steep density (n < —3) 
the structure seen in the velocity channels is determined only by ve- 
locity fluctuations (is not explicitly dependant on the index n). This 
immediately sends a warning against identifying structures in PPV 
as actual density enhancements ('clouds'). The spectral velocity in- 
dex can be obtained if we combine the density velocity index from 
the very thick channels and the structure function index m from 
thin channels (using Table Q and eq. 0. Notice that the notion of 
thin and thick slices depends on the turbulence scale under study, 
and the same slice can be thick for small-scale turbulent fluctua- 
tions and thin for large scale ones. The formal criterion for the slice 
to be thick is that the dispersion of turbulent velocities on the scale 
studied should be less than the velocity slice thickness. Otherwise 
the slice is thin. 

From TableQone may observe that for any given density spec- 
tral index, the spectrum of fluctuations becomes shallower as the 
velocity field gets steeper. This is somewhat counterintuitive. But 
can be understood as follows: a steeper velocity index corresponds 
to relative more power on the largest scales. A steeper velocity 
power-spectrum means a larger velocity difference across scales, 
and correspondingly to a larger velocity dispersion. In PPV space, 
large eddies, because their velocity dispersion, will be spread over 
a larger velocity range than smaller eddies. In individual channels, 
the number of emitters associated with large eddies will decrease 
relative to that of small eddies as the velocity contrast increases. 
Therefore, in velocity slices, the power at small k numbers (large 
scales, from large eddies) decreases for a steeper velocity, render- 
ing in a shallower power-spectrum. 

In following sections we perform tests of the VCA using nu- 
merical data, and take a step further exploring the sources of uncer- 
tainty of and the effects of shear in the VCA. 



3 TESTING THE VCA 

Numerical simulations provide good means of theory testing. Using 
numerical data, we can measure directly the velocity and density 
statistics, and we can produce synthetic spectra to analyse them 
as an observer would do. However, we need to be aware of the 



differences of synthetic and real data. For instance, only a limited 
number of points is available with numerics. This, as we show later, 
results in the most important source of uncertainty. 

3.1 The data 

Along the paper we make use of three types of data to construct the 
PPV cubes and analyse the power spectrum: 

(i) MHD simulations. We use a 3rd-order hybrid essentially- 
non-oscillatory (ENO) code to simulate fully developed turbulence. 
It was performed on a 216^ Cartesian grid, it is isothermal, com- 
pressible, and with a mean magnetic field in the x direction. To 
reduce spurious oscillations near shocks, we combine two ENO 
schemes. When variables are sufficiently smooth, we use the 3rd- 
order weighted ENO scheme (Jiang & Hu 1999) without character- 
istic mode decomposition. When oppo site is true, we use the 3rd- 
order Convex ENO scheme iLiu & Osheiiil99&) . We use a three- 
stage Runge-Kutta method for time integration. We solve the ideal 
MHD equations in a periodic box: 

|+V.(p.)=0, 

^+vVv + p-^Vidp) - (V X B) X B/47TP = /, (8) 
ot 

^ - V X (i; X i?) = 0, 

with V • S = and an isothermal equation of state. Here / is a 
random large-scale driving force, p the density, Cs the sound speed, 
V the velocity, and B the magnetic field. The rms velocity SV is 
maintained to be approximately unity, so that v can be viewed as 
the velocity measured in units of the rms velocity of the system, 
and B/y^inp as the Alfven speed in the same units. The time t 
is roughly in units of the large eddy turnover time (~ L/5V) and 
the length in units of L, the scale of the energy injection. The mag- 
netic field consists of the uniform background field and a fluctuat- 
ing field: B — Bq + b. The Alfven velocity of the mean field is 
roughly the same as the rms velocity. The average Mach number is 
~ 2.5. 

The outcomes are the density and velocity fields. The PPV cubes 
constructed with this data have a correlation between the velocity 
and the density fields which is self-consistent with MHD evolution. 
These and are used to test the applicability of the VCA where the 
basic assumption was to disregard such correlation. 

(ii) 'Modified' data. Due to numerical limitations, the power 
spectrum of the simulations mentioned above doesn't show exact 
power-law behaviour over the desired inertial range, making diffi- 
cult to test the VCA, which deals with power-law spectra. For that 
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reason we apply the same procedure as in lLazarian et alJ i200ll) to 
modify the power spectrum to follow a power-law preserving most 
of the phase correlations. The procedure consists in replacing the 
amplitudes of the Fourier transform of the data so they follow a 
power-law while keeping the phases intact (were most of the spa- 
tial information is). The set of simulations used here show a large r 
well-developed inertial range than those in lLazarian et alJ j200lh : 
therefore the 'correction' introduced is smaller, especially for the 
-11/3. 

(iii) Gaussian fields. We use a standard procedure to gen- 
erate Gaussian fields with prescribed power-law spectrum (see 
[Bond, Kofman & Pogosvan 199^. We generate a Fourier represen- 
tation of the field 

F(fc) =^fcafeP(fc)^/2e'*^ (9) 

where are independent Gaussian variables with dispersion 1, 
and a mean of for velocity fields or 1 for density fields; in our 
notation the wave-number k — 2nL / A, where A is the wavelength 
associated with k, and L is the box size. To obtain the field in con- 
figuration space a Fourier transform is applied. Velocity and density 
fields generated this way are entirely uncorrelated. However the fa- 
cility of generating as many cubes as we want in a short amount of 
time is useful when we want to improve the statistics of our sample. 

Comparing (ii) and (iii) allows to test the importance of 
density- velocity correlations. 

With the data cubes mentioned above we produced synthetic 
spectra (i. e. PPV 'cubes'^) to analyse them using VGA. The proce- 
dure to construct the PPV cubes is as follows: We generate a matrix 
of the same dimensions in X as the x data and a given number of 
velocity bins (velocity channels) as the third dimension. The ve- 
locity bins are equally spaced ranging from the minimum value of 
the velocity (vmin) to the maximum value of the velocity (vmax) 
in the whole velocity data cube. We then sum all the elements of 
the density cube in a position X that correspond to each veloc- 
ity bin. This way we simulate observational data, neglecting self 
absorption, and assuming that the thermal width of the emitting 
gas is always smaller than the thickness of the velocity channels. 
An introduction of an intrinsic width (equivalent to thermal) to the 
emitters is discussed later. 

3.2 Density and velocity correlations 

LPOO considered density and velocity either uncorrelated or cor- 
related to maximal degree. What would happen to intermediate 
cases? We computed the correlation of the density and the velocity 
fields as 

Cpau 

Where r is the spatial separation ('lag'), x is the spatial position, 
p the density field, u is the velocity field, Up and (7„ are the stan- 
dard deviations of the density and velocity fields respectively. The 
average is to be performed over all x space and angles of r. For 
numerical economy Lazarian et al. (2001) used Fourier transform 
techniques to obtain similar information (although restricted to the 
LOS velocity component). Here this is done directly using equation 
tlOt , which is applied to the original output from the MHD simu- 
lations. We found moderate correlations, the maximum correlation 

^ Strictly speaking they are parallelepipeds as the sizes need not be the 
same in all directions. 
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Figure 1. Correlation of the density and velocity fields, (a) the average pro- 
jection of the cross-correlation on the x axis (parallel to the mean magnetic 
field), (b) the average projection of the correlation on the y axis (perpendic- 
ular to the mean magnetic field), (c) the average projection of the coiTela- 
tion on the z axis (perpendicular to the mean magnetic field, and (d) average 
correlation as function of the spatial separation r. 

(in absolute value) was of ~ 0.22. This correlation is shown in Fig. 
Q where we plotted the average projection of the cross-correlation 
on the X, y, z axis and the average correlation as a function of the 
separation r. Clearly the velocity and density correlations are very 
different along and perpendicular to the magnetic field direction. 
Some implications of this anisotropy are discussed in §5. 

3.3 Applying the VGA 

In this subsection we will start the test of the VGA generating sev- 
eral data cubes with known ('modified') spectral indices, and ar- 
bitrary velocity resolutions that will be shown to range from very 
thick, to excessively thin. This is to check the general trends pre- 
dicted in LPOO. Then with a more careful choice of the velocity 
thickness of the channel, we will narrow down the measured spec- 
tral indices to check the numerical value with the analytical expec- 
tation. A more detailed analysis of the sources of error introduced 
with numerics, with an example case of Kolmogorov indices, fol- 
lows in the next subsection. 

We produced density data cubes with spectral indices modi- 
fied to —2.5 and —11/3, LOS velocity cubes with indices of —3.3, 
— 11/3 (Kolmogorov), and —4.0. With each velocity and the den- 
sity fields we constructed PPV cubes with velocity resolutions, 
ranging from 150 to 1 number of channels. The 2D power spectra 
in velocity channels (averaged over all velocity channels for each 
PPV cube) is presented in Fig.|2| 

We should stress, that whether a slice in velocity space can 
be considered thin or thick depends not only on the slice width 
5v = {vmax — Vmin)/N = Av / N (where N is the number 
of channels) but on the scale as well. For power-law statistics the 
squared dispersion, (t|,, over the box size, L, equals CI/*" (here m 
is the velocity structure function index related to the spectral index 
as described in eq|3j; similarly the velocity dispersion squared, a^, 
over the scale r is Cr™. Therefore the criteria for a channel to be 
considered thin (5v < ar) from the largest scale L up to the scale 
r translates into N > ^ {^)"^^^ ■ Thus, we do not expect a pure 
power-law result for a slice of a given width. Indeed, the largest 
scales, on the order of the size of the cube, are almost always in a 
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n = -2.5, /x=-3.3 n = -ll/3, |Ct=-3.3 
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Figure 2. Power spectra in velocity cliannels for simulated observations with density spectral indices of n = —2.5 and n = —11/3, with velocity spectral 
indices of = —3.3, 11/3, and —4 (as labelled in the title of each plot). In each panel the dotted lines con'espond to the power-spectra in velocity channels 
for PPV cubes generated with, from top to bottom 150, 100, 75, 50, 25, 15, 10, 6, 4, 2, and 1 channels respectively. For reference we plotted the predictions 
in LPOO: in dashed lines are the predictions for thin channels, in solid lines for very thick channels. 



thin regime; unless, A*' < ^ ~ 5 — 6. Figure|2lis consistent with 
this behaviour. 

In Figure |2| we can see that the power spectrum has an addi- 
tional deviation from the power-law at small scales. For that rea- 
son, to measure numerical values of the spectral indices, we have 
restricted ourselves to the wave-number range from log fc ~ 1.3 to 
2.2, where the power spectra in Fig.|2|is close to a power-law. The 
slopes (spectral indices) of the power spectra obtained with a linear 
least-squares method are shown in Table|2| 

From Fig. |2| (and Table |2j we can see, indeed, a steepening 
of the power- spectrum as we increase the thickness of the velocity 
slices. However, the spectral-indices in thin channels do not appear 
to reach an asymptotic value, instead they keep getting shallower 
as the velocity resolution increases. This behaviour arise due to the 
limited number of emitters along the LOS: with only 216 emitters 
along each LOS, to distribute in a large number of velocity chan- 
nels, there will be many channels empty. This produces artificial 
sharp edges, and introduce high frequency components. Rendering 
in a sh allower sp ectrum at large wave-numbers. To avoid this prob- 
lem, in lLazariM i et al. (2001) it was suggested as optimal to search 
for thin slice asymptotics in the channels that are just thin enough, 



at the smallest scale r of interest. The expression in eq. Jilt is 
somewhat tricky, because we need to know the value of m, to get 
the number of channels to be used to best determine m! But, it al- 
lows to check with emulated observations if the predictions hold. 
Our results, in Table |2| are in agreement with the prediction in 
[Lazarian et al. ( 2001 ) that more channels are required when steeper 
velocity fields are studied. For real observations the problem is al- 
leviated because in that case we are not limited by the number of 
emitters and we can use the maximum velocity res olution with- 
out having empty channels. In lLazarian et alj i200ll) the factor of 
Aii/cTL in eq. <1 H was estimated to be ~ 2.5, in the assumption of 
Gaussian statistics. In practice we can use the standard deviation of 
the velocity in the LOS velocity cube for ctl, and An is the range 
of velocity covered in our PPV cubes. This is important to keep in 
mind because a side effect of the modification of the spectral in- 
dex is an appreciable change in the velocity range of the velocity 
cubes. To get the factor {L/r) we need to take the maximum wave- 
number used to obtain the power-spectra. For the results presented 
in Table|2|it was k ~ 10^'^, so L/r = A:/27r ~ 25. This way 
we generated PPV cubes with the velocity resolution according to 
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Table 2. Measured emissivity spectral indices of Fig.l2l 





Density index n 


= -2.5 


Density index n = 


-11/3 


Number 






Velocity index (/x) 






of 


-3.3 


-3.67 


-4.0 


-3.3 


-3.67 


-4.0 


velocity 




Predicted spectral index (7) 




channels 


-2.35 


-2.17 


-2.00 


-2.67 


-2.6 


-2.5 


150 


-2.02 


-1.97 


-1.93 


-2.25 


-2.19 


-2.15 


100 


-2.05 


-2.00 


-1.95 


-2.27 


-2.22 


-2.18 


75 


-2.07 


-2.02 


-1.97 


-2.30 


-2.24 


-2.20 


50 


-2.11 


-2.06 


-2.02* 


-2.35 


-2.29 


-2.25 


25 


-2.22 


-2.18* 


-2.15 


-2.54 


-2.47 


-2.44* 


15 


-2.34* 


-2.33 


-2.33 


-2.74* 


-2.74* 


-2.74 


10 


-2.44 


-2.42 


-2.46 


-3.01 


-2.99 


-3.02 


6 


-2.47 


-2.52 


-2.56 


-3.13 


-3.23 


-3.39 


4 


-2.47 


-2.53 


-2.55 


-3.13 


-3.29 


-3.44 


2 


-2.46 


-2.52 


-2.54 


-3.13 


-3.30 


-3.45 


1 


-2.49 


-2.49 


-2.49 


-3.65 


-3.65 


-3.65 



* With this velocity resolution (number of channels) we get the better 
match for the analytical predictions in thin channels. 



Table 3. Emissivity spectral indices in barely thin chanels 



Density 


Velocity 


Num. of 


Predicted 


Measured 


index 


index 


channels 


7 


7 


-2.5 


-3.3 


35 


-2.35 


-2.2 


-2.5 


-3.67 


42 


-2.17 


-2.1 


-2.5 


-4.0 


54 


-2.00 


-2.0 


-3.67 


-3.3 


35 


-2.85 


-2.4 


-3.67 


-3.67 


42 


-2.67 


-2.3 


-3.67 


-4.0 


54 


-2.50 


-2.2 



eq. <llt . and measure the power spectrum in velocity channels, in 
Table |3| we display the results. From Table |3| we can see that the 
measured spectral indices within 15% error of the predictions, but 
always shallower than expected. This is an indication, that the noise 
arising from the limited number of emitters is present even in chan- 
nels barely thin, although not near as large as that present for PPV 
cubes constructed with more channels (see Table|2}. 

Could deviations above be due to correlations between the 
density and the velocity fields? The answer is no: to illustrate this 
point, we can compare the power spectrum of the MHD simulations 
to that of an ensemble of synthetic PPV cubes with the same spec- 
tral indices but no velocity-density correlations whatsoever. For this 
comparison to work we must use PPV cubes constructed from data 
with the same power-law inertial range. For that reason, we mod- 
ified the spectral index of the density field from the simulations 
to —11/3 (close to the original index, so the modification on the 
velocity-density correlations is minimal). With this density field, 
we then constructed PPV cubes using the modified velocity output 
from the simulations and 60 synthetic velocity fields. All the ve- 
locity cubes used had the same spectral index (also —11/3). The 
power spectra in velocity channels is shown in Fig. |3| In this fig- 
ure we can appreciate that apart of some 'wiggles' at small values 
of k, due to lack of enough statistics (see discussion below) the 
power-spectra is practically identical whether the density field is 
correlated with the velocity field or not. 




1.0 1.5 2.0 2.5 

log k 




-6^ (b) 

1.0 1.5 2.0 3.5 

log k 

Figure 3. Power spectrum in velocity channels for PPV produced with den- 
sity and velocity cubes from MHD simulations -with spectral indices mod- 
ified to —11/3- (a), and with the same density as the upper panel but an 
ensemble of Gaussian cubes with the same —11/3 index {b). The dotted 
lines con'espond to the average spectrum in velocity channels for PPV with 
different number of channels, in both panels from these con'espond from 
top to bottom to 150, 100, 75, 50, 25, 15, 10, 6, 4, 2, and 1 channels re- 
spectively. We see larger departures from power-law behaviour at lower k in 
(a) because in (b) we averaged over an ensemble, see discussion on the text. 
Also, in both cases for reference we plotted the predictions of the power 
spectrum in thick channels {solid lines), and in thin channels {dashed lines). 

3.4 More on the numerical limitations 

As we saw in the previous section, one important practical problem 
with the determination of the spectral indices through VGA is to 
measure with accuracy the slopes (precisely the spectral indices) in 
the plots. We can notice from each plot in Fig. |2| two departures 
from well behaved power-laws. At the left end of the spectrum 
(low k numbers, large spatial scales) we have some random wig- 
gling, while at the right end of the curves we see a systematic flat- 
tening of the power spectrum. The latest is increasingly evident for 
PPV cubes formed with larger number of channels (with thi nner ve- 
locity channels), and was referred as 'shot-noise' in iLazarian et alJ 
<200lh . 

3.4.1 Departures from power-law at low wave-numbers 

In the spirit of avoiding other sources of uncertainty, we will use 
simple Gaussian fields with prescribed power-law spectrum as ve- 
locity fields, mapped with a constant density field to produce PPV 
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cubes. We generated an ensemble of 60 Gaussian cubes with a pre- 
scribed spectral index of —11/3. Then we applied to them the 
VGA, that is, we obtained PPV cubes with different number of 
channels (different velocity resolution) and compute the 2D power 
spectra in velocity channels. The results are shown in Fig.|4| From 
the figure is evident that these departures from strict power-laws 
at low k numbers are due to the lack of sufficient statistics at the 
largest scales. Power spectra in velocity channels for individual ve- 
locity cubes differ very little from intermediate to large k numbers, 
but differ significantly at small k, showing some 'wiggles'. How- 
ever, the ensemble average converges to a power-law that extends 
from the largest spatial scales to intermediate spatial scales, until 
a systematic flattening of the power spectrum at large k becomes 
evident. This problem may ocurr in actual observations, where the 
low wave-numbers (large scales) are affected by the general size 
and shape of the object under study. 



3.4.2 Departures from power-law at high wave-numbers 

One of the most important limitations that we to construct spec- 
tral lines with numerics, is the finite (and very small) number of 
emitters in a particular direction; in our simulations we have 216 
density elements in each LOS to distribute in velocity channels. If 



we impose a high velocity resolution in the PPV cubes (i. e. large 
number of channels), the spectral lines instead of being smooth 
will have sharp edges and empty channels. This artificial jumps 
introduce high frequency components in the spectrum and there- 
fore an increase of th e power spectrum at large wave-numbers. In 
iLazarian et"ai]<200lh this was regarded as 'shot-noise'. This is true 
in the sense that the flattening becomes larger as we increase the 
velocity resolution. But is more complicated that the usual shot- 
noise, that can be easily described in terms of Poisson statistics. 
This flattening of the spectrum is particularly problematic because 
it increases as we go into thinner channels, and the predictions in 
LPOO to extract the velocity spectral index were done precisely in 
an asymptotic regime for thin channels. 

A natural way to avoid the problem of this flatening of the 
power spectrum is simply to increase the number of emitters along 
the LOS. This can be accomplished increasing the size of the data 
cubes, which translates into more computing power required to 
analyse the power-spectrum (and more importantly more comput- 
ing power to generate larger MHD cubes). With the resources at our 
hand, we increased the number of emitters along each LOS by cre- 
ating an ensemble of 50 2D Gaussian velocity fields of size 1536^. 
All with spectral index of —8/3, Kolmogorov in 2D. Analogously 
to the 3D ensemble we simulated observations and analysed the 
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Figure 5. Comparison of the power spectrum in velocity channels for an 
ensemble of 60 velocity cubes (a) with an ensemble of 50 squares (h). The 
cubes are 216'^ in size and the squares are ISSO'^. All velocity fields have 
Kolmogorov spectral indices (—8/3 for 2D and —11/3 for 3D), and the 
PPV were obtained with a constant density. The dashed line on the lower 
panel is to compare the inertial range achievable in both cases. 



power spectrum in velocity slices. In this case instead of PPV cubes 
we have simply 'PV arrays', and the power spectrum computed in 
slices is one-dimensional. In Fig.|5]we compare the ensemble av- 
erage power-spectrum in velocity channels for both, 3D and 2D 
cases. It's clear the increase in wave-number range acquired with 
the 2D data, in which we have about half decade more. The relevant 
result of this comparison is to show that the flattening of the power 
spectrum for the thinnest channels is present in both cases, at large 
k. But for the 3D case this flattening is present earlier than in the 
2D case. 

It's important to stress that this is NOT a real dynamical ef- 
fect, the cubes used to construct the PPV cubes are just Gaussian 
velocity and constant density. 

Another way to alleviate the problem of limited number of 
emitters is to include an intrinsic width to each emitting element. 
We acomplished this task by adding some Gaussian smoothing on 
each spectral line (i. e. convolving a Gaussian profile along each 
LOS). This is somewhat equivalent to introduce a thermal width 
to the emitters^. Some examples of the spectral lines in our PPV 
cubes, with and without an intrinsic width are shown in Figure|6| 



Figure 6. Examples of spectral lines in one of our PPV cubes with differ- 
ent resolutions (the number of channels is indicated in each panel.). On the 
left we have spectra from PPV cubes with no width associated to the emit- 
ting elements, on the right the same spectra after adding a significant width 
(Gaussian with FWHM of A/60). All the velocity axis contain the entire 
velocity range in the PPV, and the intensities are normalized for visual pur- 
poses. 



Introducing an intrinsic width to the PPV cubes makes the spectral 
lines to look more realistic because it smooths sharp edges due to 
the limited numerical resolution. However, it also limits the reso- 
lution that can be used to look for thin channel asymptotics. The 
criterion in LPOO for thin or thick channels, was based on an 'ef- 
fective width' of the velocity channels, which included not only the 
thickness of the channel but also a thermal width vt as: 



2n1/2 



(12) 



The additional intrinsic width limits the resolution attainable 
in two ways. For a fixed maximum k scale the velocity channels 
need to be narrower to remain in the thin regime, with a minimum 
number of channels given by: 











(i; 





-1/2 



(13) 



And for a fixed velocity resolution it reduces the wave-number up 
to which the channel remains thin to 



k < 2tt 



iV2 



(14) 



The criteria in eQ. jl4t is more restrictive, beyond that scale the ve- 
locity channels are entering in a transition to the thick regime. The 
condition in eQ. <I3> is less restrictive, in the sense that any slice 
with that minimum number of channels is to be considered thin. 
However, there is no reason to slice narrower than that, beacuse the 
thermal width smears velocity fluctuations on smaller scales. If we 
go to yet thinner channels we reduce the amplitude of the signal, 
but do not see any further change of spectral index. In Figure0we 
show the power spectra in velocity channels for PPV cubes where 
we included an intrinsic width similar to the width of the barely 
thin channels (thin down to the scale corresponding to the limit 
we chose previously to measure the power spectra, logfc ~ 2.2). 
From the figure we can notice immediately that the rising tail at 



We used this width as free parameter, therefore we rather use the term 



'intrinsic' instead of 'thermal' width, wich stricltly speaking, is determined 
by physical mechanisms beyond the scope of our simulations. 



n=-2.5, fi=-3.3 n=-ll/3, ji=-3,3 




1.0 15 S.O S5 1.0 15 i.n 'ii 

bg k bj k 



Figure 7. Same as Fig. 2 but including an intiinsic widtli to tlie emitters. 
Eacli emitting element is a Gaussian with a FWHM of Av/60. We can see 
how the flatening of the power spectrum at large k decreases dramatically. 

high wave-numbers dramatically decreases. In the cases where the 
velocity fluctuations at small scale structure are relativelly small 
(in Fig.Q n = —4) the transition from thin to thick channels as 
we go into smaller scales (large fc) is evident. This corresponds to 
predictions in LPOO. We computed the spectral indices in velocity 
channels, for our previous inertial range with slices of thickness 
according to ea. <13> . and also for the conservative inertial range 
from ea. <14t for the same velocity resolutions as before. In both 
cases we observe that the measured spectral indices doesn't change 
significantly from the earlier estimates with barely thin channels 
(presented in Table |3j- This is understandable because the intrin- 
sic width we included alleviates the problem of small number of 
emitters only on scales smaller than such width. 

This shot-like noise presents an obstacle if we want to obtain 
the velocity index of data cubes generated synthetically. However, 
we must stress that the problem is not expected to be present in 
real observations, where the number of emitters along the LOS is 
much larger. This w as confirmed through obse rvations of the Small 
Magellanic Cloud iStanimirovic & Lazarianl l?OOD. The spectral 
indices there, measured in thinner and thinner channels do tend to 
an asymptotic value, not being affected by 'shot-like noise' allowed 
to measure accurately the power-spectrum in velocity channels. 



4 EFFECT OF THE SHEAR ON THE VCA 

The VCA is sensitive to density and velocity fluctuations at a partic- 
ular scale. On the other hand, shear as any organized motion, like 
galactic rotation, would correspond to the smaller wave-numbers 
(the largest scales, sometimes much larger that the largest scale in a 
PPV cube). If we choose the right range in k to perform the VCA, 
we should avoid any interference of this ordered motions with the 
analysis. To explore numerically the effect of a linear shear on the 
VCA we used modified data cubes of spectral indices of —11/3 for 
both velocity and density fields. We chose the LOS to be along the 
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Table 4. Spectral indices in velocity channels with shear 



Num. of No Shear magnitude (in units of (T„^ /L) 
channels shear 1 3 5 10 20 



150 


-2.2 


-2.3 


-2.4 


-2.4 


-2.7 


-3.1 


100 


-2.3 


-2.3 


-2.4 


-2.5 


-2.7 


-3.1 


75 


-2.3 


-2.4 


-2.5 


-2.5 


-2.8 


-3.2 


50 


-2.4 


-2.4 


-2.6 


-2.6 


-2.9 


-3.3 


25 


-2.6 


-2.7 


-2.8 


-2.9 


-3.2 


-3.7 


15 


-2.9 


-3.0 


-3.0 


-3.1 


-3.4 


-4.1 


10 


-3.2 


-3.1 


-3.1 


-3.3 


-3.4 


-4.0 


6 


-3.4 


-3.3 


-3.1 


-3.1 


-3.3 


-3.7 


4 


-3.4 


-3.3 


-3.1 


-3.1 


-3.3 


-3.5 


2 


-3.4 


-3.3 


-3.1 


-3.1 


-3.2 


-3.4 


1 


-3.7 


-3.7 


-3.7 


-3.7 


-3.7 


-3.7 



z axis and added a component to the velocity cube (the cube) 
following 

X 

Vz{x,y,z) = voz{x,y,z) + Cav^yZ, (15) 

where is the new velocity cube used to generate the PPV cubes, 
voz is the original velocity cube before introducing the shear, a^^ is 
the velocity dispersion of voz , L is the cube size in the spatial coor- 
dinates (x, y), and C is a parameter we introduce to vary the mag- 
nitude of the shear introduced. The shear introduced is the same on 
all scales, and equal in magnitude, to Ca-v^ /L. We generated PPV 
cubes adding a linear shear with magnitude up to 1, 3, 5, 10, and 
20 times maximum velocity dispersion over the largest scale (i. e. 
with a velocity component with C = 1, 3, 5, 10, and 20 in eq. 
I15t . As before, we measure the spectral indices in the range less 
affected by numerical noise, log k from ~ 1.3 to ~ 2.2. The range 
used is somewhat arbitrary, but the purpose of the study at this point 
is just comparison with the case of no shear The measured slopes 
are summarized in Table|4| The introduced shear is the same on all 
scales, while the turbulent shear scales as ~ k^^^^^^^K The mag- 
nitude of the external shear is the same as the magnitude of the tur- 
bulent shear when k ~ 27rC"^''^. Therefore the external shear con- 
tributed initially to the power spectrum only at the largest scales, 
but eventually became important at intermediate to small scales as 
its magnitude increased (see Fig.|8}. Our results show that the effect 
of an external shear becomes important only when it is comparable 
to the local turbulent shear. In our example that was when the max- 
imum shear reached about 5 times the velocity dispersion over the 
largest scale. The relative differences of the spectral indices mea- 
sured compared with the case of no shear are below 10%, for the 
cases of less shear (C < 5). It is noticeable that the larger distor- 
tion happens at lower k, even when its magnitude is comparable to 
the turbulent shear at smaller scales. This is expected, given that the 
shear is regular, and represents a large scale motion. And holds if 
we measure the slope of the power spectrum in intermediate spatial 
scales, which are available using synthetic data. With actual obser- 
vations, we can go further into smaller scales and avoid the effects 
of a stronger shear. The shear due to turbulence dominates in the 
galaxy, and its the effect to the VCA is marginal. 



5 ANISOTROPIC TURBULENT CASCADE 

Magnetic field plays a crucial role in the dynamics of inter- 
stellar turbulence, and makes the turbulent cascade anisotropic 
( lMontgomer^ll98llHigd oi]|l984) . In a turbulent magnetized me- 
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uding a shear of 1 x o^^L 




Figure 8. Power spectra in velocity cliannels including linear shear Each 
panel correspond to a different magnitude of the maximum shear, as labelled 
in the title of each plot. The different lines in each panel correspond from 
top to bottom to the average power spectrum in velocity channels for PPV 
cubes constructed with 150, 100, 75, 50, 25, 15, 10, 6, 4, 2, and 1 channels 
respectively. The small aiTows show at which wave-number the magnitude 
of the external shear equals the magnitude of the turbulent shear. 



dia, the energy on the large scale motions is larger than that on the 
small scale structures; however the local magnetic field strength 
remains almost the same. Thus, is easier to bent magnetic field 
lines on the large scale by turbulent motions, but on the small 
scale becomes more difficult. Another way to think about this 
phenomenon is noting that hydrodynamic motions can easily mix 
magnetic field lines on the small scale, rather than bent them, in 
the directi on perpendicular to t he mean magnetic field (see a dis- 
cussion in iLazarian & VishniaQ. 199 9 V As a result we get elon- 
gated eddies, relative to the magnetic field that become even more 
elongated as we go further into smaller scales. A model for in- 
compr essible anisotropic turbulence is that ofjjgldreich & Sri dhar. 
Jl99a). This model has been supported by numerical simu- 



lations (e.g. |C ho & Vis hniad [20001 : 'Maron & GoldreiclJ 1200 it 
ICho. Lazarian & Vishniacl l2002 ). Compressible anisotropic MHD 



Lithwick & GoldreicM i2001h and 



turbulence has been studied in 
Icho & Lazariaa(.2002,.2003l) . 



5.1 Magnetic Field Direction: Anisotropy in Correlation 
Functions 

In an isotropic turbulent cascade correlations depend only on the 
separation between points. Contours of equal correlation are circu- 
lar in that case. If on the other hand the turbulence is anisotropic, 
two point statistics are also anisotropic. In this case contours of 
iso-correlation become elongated, with symmetry axis given by the 
magnetic field direction. The technique presented here is similar to 




5 10 15 20 25 50 

Distance parallel to ttie B tield 



Figure 9. Contours equal correlation of 2D maps obtained with the PPV 
data, (a): contour map of the correlation of emissivity at the centroid of 
the lines, (h): contours of correlation in the integrated emissivity. (c), (d): 
contours in two individual channel maps. These were obtained in a PPV 
cube with 25 channels (the velocity range is from Vmin = —1.27 to 
Vmax = 1.451, and consequently the channel thickness is SV = 0.109). 
(c) corresponds to a velocity centred at Vz = —1.11, while the (d) to a 
velocity of Vz = 0.741. All the distances in the figures are in grid units. 



one proposed to study magnetic field direction using synchrotron 
maps ILazarian 1992). Iso-contours of correlation functions can al- 
low us to determine the direction of the magnetic field. We calculate 
two point correlation for different 2D maps generated with our data. 
The correlation function between the two scalar functions f{x) and 
g{x) can be computed from 

{{f{x)-{f{x))).{g{x + r)-{g{x)})} 



(TjCJg 



(16) 



where r is the separation between two points or 'lag', x is the spa- 
tial separation, a/ and (Jg are the standard deviations of / and g 
respectively. Here the average is performed on all the radial com- 
ponents and azimuthal angles. 

With the original cubes from simulations (no artificial modi- 
fication of the spectral indices), we calculated the auto-correlation 
(correlation in eq. jl6> with f{X) — g{X)) on various 2D maps: 
integrated intensity (integrated along the LOS, equivalent to a 2D 
map of column density); emissivity at the centroid of the lines, ob- 
tained as: 



MX) 



Jj{X)vdv 



(17) 



JjiX)dv ' 

where the emissivity j{X) in our case is proportional to the first 
power of the density; and maps of emissivity in individual veloc- 
ity channels. Contour maps of iso-correlations calculated of these 
maps are given in Fig.|9| We can see that anisotropy is present in 
the contours plotted. It is evident in the case of the emissivity at the 
centroids of the lines, as well as in the integrated emissivity. In in- 
dividual velocity slices, we found that the majority of the velocity 
channels show this anisotropy in the expected direction (elongation 
of the contours parallel to the mean B field), but in some slices 
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we see the elongation of correlation contours tilted from the mean 
magnetic field direction. For Fig.|9|we choose a representative slice 
to illustrate this misleading behaviour. 

The main purpose of this subsection was to show an alternative 
method to obtain the direction of the magnetic field, and also to 
illustrate that anisotropic cascade is present in our simulations. 

The anisotropy based technique could be a unique tool for 
magnetic field studies, especially when other techniques, e.g. based 
on dust alignment fail (see iLazariaa.200Qi for a review on grain 
alignment). 

5.2 Anisotropic cascade & VCA. 

It's important to test the VCA, which was developed considering 
isotropic turbulence, in a more realistic case where the magnetic 
field determines a preferential direction of the turbulent motions. 
With that purpose we constructed PPV cubes choosing the LOS to 
be parallel to the mean magnetic field (along the x axis) and com- 
pare the results with the case where the magnetic field is perpendic- 
ular to the LOS (z axis). For this exercise we use a data cubes cor- 
rected to Kolmogorov spectral indices (—11/3 in 3D). The power 
spectrum in velocity channels is shown in Fig. llOl In the figure we 
can see that the slopes of the spectrum (spectral index) are almost 
identical whether the observer is looking parallel or perpendicular 
to the mean magnetic field. We performed a least-squares linear fit- 
ting and found the differences of the slopes in every case are well 
within the fitting uncertainties. 

We conclude that VCA is applicable to anisotropic turbulence 
studies. 



6 SUMMARY AND CONCLUSIONS 

In this paper we have tested the VCA technique numerically. We 
used compressible MHD simulations, and an ensemble of Gaus- 
sian fields to emulate spectroscopic observations (PPV data cubes). 
We computed the power spectra in velocity channels to test the an- 
alytical predictions in LPOO. 

To obtain the PPV data we modified the output of the sim- 
ulations to force power-law indices to obtain density fields with 
spectral indices of —2.5 and —11/3, and velocity fields with spec- 
tral indices of —3.3 to —4.0). The spectral indices derived from our 
emulated observations are within a 15% error of the analytical pre- 
dictions. We confirmed trends predicted in LPOO, like the steepen- 
ing of the emissivity index as we increase the width of the velocity 
chaimels, or that the velocity fluctuations dominate the emissivity 
for thin velocity channels. This last result is in agreement with a 
previous study by Pichardo et al. (2000), and is important because 
warns observers against interpret any structure in PPV as density 
enhancements or 'clouds' (see also discussion in LPOO). 

We showed that the deviations from power-law behaviour at 
the largest spatial scales (small k numbers) are due to lack of good 
statistics at those scales. This was comparing the results obtained 
with PPV produced with a constant density and the modified veloc- 
ity field of spectral index — 1 1/3, to those using an ensemble of 60 
Gaussian cubes with the same index. 

An important source of error addressed is a flattening of the 
power spectrum at large wave-numbers. This is due to the lim- 
ited number of emitters along the LOS. We illustrated this point 
comparing the results obtained with the ensemble of 60 cubes with 
216 emitters along the LOS, and an ensemble (two-deimensional) 
with 1536 emitters. We also try including an intrinsic width to the 



emitters along the LOS to smooth the sharp edges arising from 
the limited number of emitters. This procedure reduced the flat- 
tening of the spectrum at large k numbers. But, also increases the 
effective width of the channels, reducing our ability to make use of 
the scales on wich we reduced such flattening of the power spec- 
trum. This problem was found to be the most restrictive to our 
work. However, is exclusive of numerical simulations, and is not 
expected to be important in real observations (like those presented 
in Stanimirovic & Lazarian 2001), where the number of emitting 
elements is essentially infinite. 

We introduced linear shear to our data, and showed that its 
effect only becomes important when the magnitude of the shear is 
comparable to the local turbulent shear. We showed, that the effect 
of large scale shear is marginal at small scales. And that avoid- 
ing the largest scales to computed the power-spectrum allowed to 
recover the turbulence statistics. This shows that for ordinary ob- 
servational situations we can safely neglect shear in the VCA. 

We analysed the cross-correlation for 2D maps of integrated 
intensity and emissivity at the centroids of the lines, and showed 
that they are anisotropic (elongated), with symmetry axis defined 
in the direction of the mean magnetic field. This constitutes a 
new technique to study the direction of the magnetic field. We 
hope that the anisotropics will reveal the magnetic field within 
dark clouds, where grain alignment and therefore polarimetry fails. 
From the point of view of VCA, we illustrated that the effect of the 
anisotropy of the turbulence is marginal. 
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